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B  tjjg  operation  and  use  of  a  new  numerical  model, 

GNWave.  This  model  was  designed  to  simulate  the  evolution  of  a  train  of  two- 
dimensional  waves  in  waters  of  arbitrary  bottom  topography,  varying  from  shallow 
water  to  waters  of  moderate  depth.  The  program  uses  the  Green-Naghdi  theory  of 
fluid  sheets  as  Its  model,  and  Integrates  a  set  of  coupled,  nonlinear  partial 
differential  equations  in  time  to  perform  the  simulation  of  surface  gravity 
waves .  Report  1  In  this  series ,  entitled  "Application  of  the  Green-Naghdi  theory 
of  fluid  sheets  to  shallow  water  wave  problems,"  contains  a  detailed  description 
of  the  mathematical  basis  of  GNWave  model.  The  model  GNWave  has  been  shown  to 
reproduce  with  engineering  accuracy  the  evolution  of  a  wave  of  permanent  form, 
from  small  amplitudes  up  to  almost  breaking  conditions.  The  numerical  model  Is 
portable  with  little  or  no  changes  to  a  wide  variety  of  platforms,  and  has 
successfully  been  tested  on  high-end  PC's  and  VAX  and  CRAY  mainframe  systems. 

The  governing  equations  were  programmed  using  Fortran  as  the  language.  The 
program  consists  of  a  main  program  and  a  number  of  modules  to  perform  the 
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calculations.  A  functional  flow  chart  of  individual  routines  involved  in  the 
computation  is  included  herein.  The  main  routine  directs  the  sequence  of  the 
calculation,  Including  the  reading  of  the  input  and  integration  of  the  equations, 
and  performs  some  post-processing.  The  main  program  is  divided  into  two  basic 
segments :  an  input  portion  where  the  specifications  of  the  problem  are  read  in 
and  echoed  to  the  computer  screen,  and  an  integration  portion  where  the  equations 
are  Integrated  one  time  step  at  a  time.  Output  is  produced  at  pre- specif led 
instants  in  the  computation  stages. 

The  calculation  involves  four  basic  components:  generation  of  waves  at  the 
ocean  end  of  the  wave  tank,  enforcement  of  an  appropriate  boundary  condition  at 
the  shore  end  of  the  computational  domain,  solution  of  the  two -point  boundary 
value  problem  in  space  at  a  given  instant  in  time ,  and  integration  of  the 
equations  in  time.  Each  of  these  four  components  requires  special  techniques, 
which  are  discussed  in  this  report  to  give  the  user  an  overview  of  the  flow  of 
information  and  a  flavor  for  the  computation  that  occurs  during  a  simulation. 
To  further  assist  the  reader  in  understanding  the  Implementation  of  the  code,  the 
algorithms  used  for  wave  generation,  the  treatment  of  beach-end  boundary 
conditions ,  and  the  spatial  and  temporal  integration  of  the  governing  equations 
are  explained  in  the  body  of  this  report.  Model  input  and  output  features  are 
also  presented,  together  with  the  procedure  for  using  the  program  and  the  re¬ 
start  mode  of  operation  of  the  model. 

Five  example  problems  typically  encountered  in  the  military  and  civil  works 
areas  of  the  US  Army  Corps  of  Engineers  are  presented.  These  include  collision 
of  waves  with  a  hydraulic  spillway  gate,  waves  transforming  over  an  open-ended 
sloping  beach,  pseudo-spectral  waves  passing  over  a  sand  bar,  random  waves 
propagating  through  a  dredged  navigation  channel,  and  regular  waves  shoaling  on 
a  landing  beach  of  uniform  slope.  The  first  example  problem  tests  the  model's 
capability  for  a  reflective  boundary  while  the  second  problem  checks  the 
Implementation  of  the  transmissive  or  open  type  boundary  condition.  Simulation 
results  for  problems  3  and  4  are  compared  to  illustrate  the  effects  of  a  sand  bar 
on  the  down-wave  progressing  in  a  channel.  The  last  problem  provides  an 
opportunity  to  check  the  model  simulation  against  physical  model  test  data  from 
a  laboratory  study  conducted  at  the  Danish  Hydraulic  Institute.  The  model  has 
been  shown  to  reproduce  with  engineering  accuracy  the  evolution  of  regular  and 
irregular  wave  trains  up  to  the  almost  breaking  condition. 
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PREFACE 


This  research  study  was  conducted  jointly  by  the  Coastal  Engineering  Research  Center 
(CERC),  US  Army  Corps  of  Engineers  Waterways  Experiment  Station  (WES),  and  Ship 
Research  Incorporated  (SW),  affiliated  with  the  University  of  California  at  Berkeley.  The  study 
was  authorized  under  the  Time-Dependent  Nonlinear  Free-Surface  Wave  Modeling  Project  as 
part  of  the  In-House  Laboratory  Independent  Research  Program  (ILIR),  Project  A91D,  in 
accordance  with  guidelines  established  by  the  WES  Office  of  Technical  Planning  and  Programs. 

The  study  was  performed  and  the  report  prepared  over  the  period  1  May  1991  through  1 
September  1992  by  Dr.  Zeld  Demirbilek,  Research  Division  (RD),  Coastal  Oceanography 
Branch  (COB),  CERC,  and  Dr.  William  C.  Webster,  Professor  of  Naval  Architecture  and  Ocean 
Engineering  at  the  University  of  California,  Berkeley  and  President  of  Ship  Research 
Incorporated  of  Kensington,  California.  The  work  was  conducted  under  Purchase  Order 
DACA39-91-M-5515.  This  is  the  second  in  a  series  of  reports  presenting  the  operation  and  use 
of  a  new  numerical  model  to  simulate  the  evolution  of  surface  gravity  waves. 

This  study  was  performed  under  the  general  supervision  of  Dr.  James  R.  Houston  and  Mr. 
Charles  C.  Calhoun  Jr.,  Director  and  Assistant  Director,  CERC,  respectively;  Mr.  H.  Lee 
Butler,  Chief,  RD,  CERC;  and  under  the  direct  supervision  of  Dr.  Martin  C.  Miller,  Chief, 
COB,  CERC.  WES  Project  Manager  during  the  conduct  of  this  study  and  the  publication  of  this 
report  was  Ms.  Mary  Vincent. 

At  the  time  of  publication  of  this  report,  Director  of  WES  was  Dr.  Robert  W.  Whalin. 
Commander  and  Deputy  Director  was  COL  Leonard  G.  Hassell,  EN. 
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This  report  is  aimed  at  describing  the  operation  and  use  of  a  new  computer  program, 
GNWave.  This  program  was  designed  to  simulate  the  evolution  of  a  train  of  two-dimensional 
waves  in  waters  of  arbitrary  bottom  topography,  varying  from  shallow  water  to  waters  of 
moderate  depth.  The  program  uses  the  Green-Naghdi  theory  of  fluid  sheets  as  its  model,  and 
integrates  these  coupled,  nonlinear  partial  differential  equations  in  time  to  perform  the 
simulation.  The  model  has  been  shown  to  reproduce  with  engineering  accuracy  the  evolution 
of  a  wave  of  permanent  form,  from  small  amplitudes  up  to  almost  breaking  conditions.  The 
program  is  written  in  the  Fortran  language  and  is  portable  with  little  or  no  modification,  to  a 
wide  variety  of  platforms.  The  program  development  work  was  performed  on  an  Apple 
Macintosh  nfx  and  the  program  can  run  successfully  on  other  high-end  PC’s  or  on  mainframes. 


Figure  1.  A  typical  flow  for  simulation  by  GNWave 

Figure  1  illustrates  a  typical  physical  situation  which  can  be  investigated.  The 
program  is  therefore  a  nonlinear  numerical  wave  tank  in  which  the  bottom  profile  can  be 
arbitrarily  chosen  and  in  which  wave  gauges  can  be  positioned  at  will.  The  types  of  studies 
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which  can  be  performed  are  varied  and  include  problems  of  both  military  and  civil  interest. 
One  can  study,  for  instance,  the  effect  of  a  submerged  sandbar  on  a  train  of  waves 
approaching  a  beach  (of  critical  importance  during  military  landings),  the  reflection  of  waves 
and  the  forces  on  a  seawall  or  sluice  gate  (important  for  the  design  of  civil  works  projects), 
or  the  time  history  of  a  wave  pressure  gage  under  a  train  of  near-breaking  waves  (of  interest 
to  scientific  researchers). 

A  great  deal  of  care  was  taken  in  the  programming  of  this  algebraically  intensive 
theory.  The  code  was  written  using  the  very  rigorous  scientific  programming  standards 
established  by  Ship  Research,  Incorporated,  and  comments  and  documentation  are  included  in 
the  code  itself.  The  code  was  given  extensive  tests  and  is  believed  to  be  bug-free,  but  with  a 
code  of  this  complexity,  it  is  not  feasible  to  exercise  every  possible  logic  path  and,  therefore, 
the  possibility  exists  that  some  well-hidden  bugs  remain.  The  intention  was  to  develop  a 
robust  code  which  can  be  used  widely.  Care  was  taken  to  make  the  input  and  output  as 
simple  and  as  useful  as  possible. 

The  sections  below  outline  the  GNWAVE  code  itself,  including  a  brief  description  of 
the  Green-Naghdi  theory  which  forms  the  basis  of  this  program,  as  well  as  input  and  output. 
Several  worked  examples  are  presented  to  demonstrate  the  uses  of  the  program  and  its 
versatility  for  practic^  problems  of  interest  to  the  Corps. 
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The  theoretical  basis  for  this  computer  program  are  the  Green-Naghdi  Level  n 
shallow-water  equations.  These  equations  were  originally  constructed  by  Shields  (1986)  and 
reported  by  Shields  and  Webster  (1988).  In  these  original  works,  the  results  of  the  theory 
were  successfully  compared  with  large-amplitude  waves  of  permanent  form  and  with  some 
experiments  on  wave  shoaling  (these  latter  results  are  reproduced  here  as  a  worked  example). 
In  the  original  developments,  the  equations  were  expressed  with  a  vertical  coordinate  which 
was  non-dimensional.  In  the  study,  these  equations  were  re-derived  in  dimensional  form, 
making  them  more  amenable  for  use  in  the  general  simulation  code  developed  here 
(Demirbilek  and  Webster  1992). 

The  Green-Naghdi  theory  of  fluid  sheets,  hereafter  referred  to  simply  as  GN  theory, 
was  first  proposed  over  20  years  ago.  It  is  an  approach  which  is  fundamentally  different 
from  the  developments  in  classical  wave  theory  begun  by  Stokes  and  Boussinesq  in  the  last 
century.  In  classical  theory,  a  perturbation  expansion  is  developed  in  which  the  nonlinear 
boundary  conditions  are  simplified  using  an  assumption  that  some  characteristic  of  the  wave 
is  small  (usually  the  wave  amplitude).  The  field  equation  (Laplace’s  equation)  is  linear  and 
does  not  require  such  an  expansion.  For  waves  of  permanent  form  this  results  in  a  sequence 
of  linear  problems  of  increasing  complexity,  which  can  be  solved  to  determine  the  flow. 

Since  each  problem  can  be  solved  separately,  higher-order  solutions  (such  as  Stokes  5*  order 
solutions)  are  commonly  used.  For  unsteady  problems,  such  as  those  treated  here, 
perturbation  expansions  yield  evolution  equations  that  become  complex  and  nonlinear  with 
each  new  order.  Few  researchers  have  performed  computations  with  equations  of  order 
higher  than  the  lowest  order  that  has  nonlinear  terms. 

GN  theory  is  based  on  an  approach  different  from  the  classical  approach.  In  this 
approach,  the  form  of  the  velocity  variation  in  the  vertical  direction  across  the  fluid  sheet  is 
assumed  to  be  given.  The  horizontal  coordinate  is  taken  to  be  x  and  the  vertical  coordinate 
is  z.  In  the  Level  n  equations  used  in  the  program,  the  horizontal  velocity  is  assumed  to 
vary  linearly  in  z  and  the  vertical  velocity  quadratically  in  z.  Thus,  for  any  x,  it  is  assumed 
that  the  horizontal  velocity  is  given  by 

u(x,z,t)  =  u.(x,t)  +  u,(x,0  z 

whereas  the  vertical  velocity  is  expressed  as 

v(x,z,t)  »  v.(x,t)  +  v,(x,t)  z  +  v,(x,t)  z^ 

The  equations  of  motion  or  governing  equations  are  derived  by  enforcing  exact 
kinematic  and  dynamic  boundary  conditions  on  the  free  surface  and  on  the  bottom,  as  well  as 


conservation  of  mass  and  conservation  of  the  O'*  and  1“  moments  of  momentum  across  the 
sheet  (in  the  vertical  direction).  These  conditions  yield  eleven  coupled  partial  differential 
equations,  which  can  be  reduced  to  three  rather  complicated  governing  equations  by 
elimination  of  many  of  the  variables.  The  remaining  dynamic  variables  are:  ;?(x,t),  u„(x,t) 
and  u,(x,t).  These  final  equations  and  their  derivation  can  be  found  in  Demirbilek  and 
Webster  (1992). 

In  summary,  the  GN  theory  is  different  from  the  perturbation  approach  in  that  the 
free  surface  and  bottom  boundary  conditions  are  met  exactly,  whereas  the  field  equation  is 
implicitly  approximated.  The  result  is  a  theory  that  can  pr^ict  the  shape  and  behavior  of 
waves  up  to  almost  breaking  conditions,  at  least  for  those  cases  of  waves  of  permanent  form 
where  numerically  accurate  solutions  are  known.  GN  theory  breaks  down  when  the  particle 
velocity  at  the  crest  equals  the  wave  speed,  the  criterion  for  breaking  in  the  exact  theory. 
Limited  comparisons  with  shoaling  wave  tests  have  also  shown  that  this  theory  can  reproduce 
this  complicated  process. 


Computer  Code 

The  governing  equations  were  programmed  using  Fortran  as  the  language.  The 
program  consists  of  a  main  program  and  a  number  of  subroutines  to  perform  the  calculation. 
A  functional  flow  chart  of  the  individual  routines  involved  in  the  computation  is  shown  in 
Figure  2  on  the  following  page.  The  main  routine,  GNWave ,  directs  the  sequence  of  the 
calculation  including  the  reading  of  the  input,  integration  of  the  equations,  and  performing 
the  output.  The  GNWave  program  is  divided  into  two  basic  segments:  an  input  portion 
where  the  specifications  of  the  problem  are  read  in  and  echoed  to  the  computer  screen,  and 
an  integration  portion  where  the  equations  are  integrated  one  time  step  at  a  time.  Output  is 
produced  in  this  loop  at  various,  pre-specified  instants  in  the  computation.  The  remainder  of 
this  section  focuses  on  the  computational  aspects  of  GNWave;  the  input  and  output  are 
discussed  in  separate  sections  below. 
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GNWave 
(Main  Program) 


Subroutines 
called  from 
GNWave 


AuxillarySubroutines 
&  functions  called  by 
other  subroutines 


Figure  2.  Program  level  flow  chart  for  GNWAVE 


The  calculation  involves  four  basic  components:  generation  of  waves  at  the  ocean  end 
of  the  wave  tank,  enforcement  of  an  appropriate  boundary  condition  at  the  shore  end  of  the 
tank,  solution  of  the  two-point  boundary  v^ue  problem  in  space  (at  a  given  instant  in  time), 
and  integration  of  the  equations  in  time.  Each  of  these  four  components  requires  special 
techniques,  which  are  discussed  below.  The  purpose  for  this  discussion  is  to  give  the  user  an 
overview  of  the  flow  of  information  and  a  flavor  for  the  computation  that  occurs  during  the 
simulation.  It  is  not  intended  to  be  a  step-by-step  presentation  of  the  algorithm.  The 
otation  used  below  corresponds,  wherever  possible,  to  that  used  in  the  actual  Fortran  code. 

To  assist  the  reader  in  understanding  the  implementation  of  the  code,  the  algorithms 
sections  of  the  code  associated  with  wave  generation,  the  right-hand  boundary  condition,  the 
spatial  integration  of  the  governing  equations,  and  the  temporal  integration  of  the  governing 
equations  will  be  described  next. 
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PARTm:  MODEL  ALGORITHMS 
Wave  Generation 


At  the  beginning  of  a  simulation,  the  water  in  the  numerical  wave  tank  is  quiescent 
and  there  are  no  waves.  It  is  noted  in  passing  that  this  is  the  only  known  exact  solution  for 
nonlinear  waves  passing  over  an  uneven  bottom.  The  computational  domain  is  discretized 
into  equal  spatial  steps  of  length  dx,  leading  to  ns  points.  The  dx  is  selected  so  that  at  the 
water  depth  corresponding  to  the  origin  of  the  computational  domain  each  wave  length  is  100 
dx  if  regular  waves  are  specified;  or  the  shortest  wave  is  60  dx  if  irregular  waves  are 
speciOed.  At  each  point  in  the  computational  domain,  spatial  derivatives,  up  to  the  third 
derivative,  need  to  be  determined  for  the  independent  variables'.  To  accommodate  this 
need,  two  additional  points  are  added  to  the  left  of  the  origin  of  the  computational  domain 
and  two  more  points  are  added  to  the  right  end  of  the  computational  domain.  The  points  are 
indexed  from  i  =  1  to  ns  +  4,  where  the  computational  domain  corresponds  to  the  points 
from  i  =  3  (at  X  =  0)  to  i  =  ns  +  2  (at  X  =  1)^. 

Waves  are  to  be  imposed  on  the  left-hand  boundary  of  the  domain  and  subsequently 
propagate  to  the  right.  In  order  to  facilitate  this  process,  the  following  approach  is  adopted. 
The  time  discretization,  dt,  is  selected  so  that  c  dt  =  dx,  where  c  is  the  regular  wave 
celerity  (for  a  water  depth  the  same  as  at  the  wave  maker)  or  the  celerity  of  the  wave  of 
median  frequency,  if  an  irregular  wave  train  is  selected.  With  these  selections,  time  histories 
for  each  of  the  three  independent  variables  are  constructed.  In  the  case  of  regular  waves,  the 
subroutine  mkwav  is  called.  This  subroutine  calls  two  other  subroutines,  wavubn  and  wvspbc. 
WAVLEN  estimates  the  wave  length  of  a  regular  wave  of  this  height  and  period  in  the  water 
depth  at  the  wave  maker;  wvspec  interpolates  a  series  of  tables  to  determine  the  first  five 
Fourier  coefficients  of  P,  u,,  and  u,  for  regular  waves  of  this  height  and  length. 

These  data  are  used  by  mkwav  to  construct  the  required  time  histories  of  the  three 


'  The  bottom  profile,  o(x),  must  also  have  third  derivatives.  As 
described  later  in  this  section,  the  subroutine  BOTTOM  constructs  such  a 
bottom  profile. 


^  The  actual  wave  maker  boundary  condition  is  applied  to  the  point  i  >  n, 
+  1  and,  thus,  the  point  i  =  n,  +  2  is  the  first  point  in  the  fluid  field. 
Similarly,  the  open  boundary  condition  is  applied  at  i  =  n,  +3,  with  the 
point  i  =  n,  +  2  being  the  last  point  in  the  fluid  field.  In  the  case  of  a 
reflection  boundary,  the  point  i  =  n,  +  2  is  the  point  of  reflection  and 
this  is  a  legitimate  fluid  field  point. 
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independent  variables  for  t  =  0,  dt,  2  dt ...  nt  dt,  where  nt  dt  =  tmax.  Because  the  fluid 
field  is  initially  quiescent,  the  wave  time  history  is  “ramped”  so  that  the  wave  maker 
imposes  only  a  little  disturbance  at  first  and  gradually  builds  to  the  full  disturbance.  The 
ramp  is  accomplished  by  multiplying  the  time  history  from  t  =  0  to  t  =  T  (the  period  of  the 
regular  wave  or  the  period  of  the  median  frequency  wave  of  a  spectra)  by  a  bell-shaped 
function  formed  by 

sinH/r  t  /  (2  T)] 

This  function  is  0  and  its  slope  is  0  at  t  =  0,  and  is  1  and  its  slope  is  0  when  t  =  T.  This 
ramp  function  forms  time  histories  which  are  continuous  and  smooth,  even  at  t  =  T. 

By  construction,  this  wave  time  history  is  identical  with  the  desired  spatial  history  for 
regular  waves  and  is  nearly  identical  for  the  case  of  irregular  waves,  especially  if  the 
spectrum  is  relatively  narrow-banded.  For  example,  for  waves  of  permanent  form  moving  in 
the  positive  x-direction,  we  have 

=  0(\-c  t) 

and  with  the  choice  of  dx  =  c  dt,  we  then  have 
^[(x„  -H  dx)  -  c  tj  =  )S[x<,  -  c  (t„  -  dt)] 

The  imposed  time  history  can  therefore  be  used  to  define  ;ff(x,t),  u„(x,t)  and  u,(x,t)  at  i 
=  1,2  and  3  (corresponding  to  the  left-hand  edge  of  the  computational  domain  and  the  two 
points  outside  of  the  domain).  These  points  provide  a  boundary  value  for  the  GN  equations 
at  this  end  of  the  computational  domain.  At  each  time  step,  the  values  of  the  three 
independent  variables  at  i  =  1,2  and  3  are  updated.  This  process  simulates  a  wave  maker 
on  the  left-hand  edge  of  the  computational  domain. 


Right-hand  Boundary  Condition 

GNWave  is  programmed  to  simulate  two  different  boundary  conditions  on  the 
right-hand  side  of  the  domain:  a  reflection  boundary  condition  and  an  open  boundary 
condition.  If  the  waves  are  perfectly  reflected  from  the  right  hand  boundary,  then  we  require 
that  the  boundary  at  i  =  ns  +  2  act  as  a  fluid  “mirror”.  For  instance,  we  set 

fi\m  +9  ~  P\m*  1 

for  the  points  outside  of  the  fluid  domain  on  the  right  hand  side,  and  likewise 
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However,  similar  conditions  must  also  be  set  on  the  horizontal  velocity,  that  is 


U.L. 


and  u.L^, 


and 


u.L  +  3  =  -u,|,  +  „and  u,|,*4  =  -uj_. 

That  is,  at  the  reflection  boundary  the  wave  amplitude  is  symmetric  and  the  fluid  velocity  is 
anti-symmetric  (resulting  in  zero  fluid  velocity  at  the  reflected  boundary).  The  reflection 
boundary  as  stated  above  is  exact  if  the  boundary  is  vertical  and  applies  equally  to  linear  and 
nonlinear  waves. 

The  simulation  of  an  open  boundary  (that  is,  the  simulation  of  the  situation  v/here  the 
waves  continue  to  propagate  unfettered  beyond  the  computational  domain)  is  much  more 
difficult.  No  exact  equivalent  for  this  boundary  condition  is  known  for  nonlinear  water 
waves.  In  other  simulations,  a  variety  of  different  methods  have  been  used.  For  the 
puiposes  of  this  study,  the  Schrodinger  condition  is  adopted,  where  the  wave  celerity  is 
assumed  known  and  constant  at  the  right-hand  boundary.  In  this  case,  one  can  write  simple 
expressions  for  the  values  of  the  three  independent  variables  at  the  right-hand  edge  of  the 
computational  domain  and  the  two  points  just  outside  of  the  domain. 

Experience  has  shown  that  this  condition  simulates  an  open  boundary  condition  for  a 
considerable  interval  of  time  after  the  wave  train  impinges  on  this  boundary.  However, 
experience  also  shows  that  if  the  waves  are  nonlinear,  this  condition  eventually  breaks  down 
and  leads  to  instabilities,  lack  of  conservation  of  mass  in  the  computational  domain,  or 
anomalous  wave  breaking.  There  appears  to  be  no  rule  of  thumb  which  predicts  when  this 
breakdown  will  occur.  In  forming  the  open  boundary  condition  for  GNWave,  the  nonlinear 
celerity  of  waves  of  the  same  period  and  height  as  the  incoming  waves  is  computed  for  the 
water  depth  where  the  condition  is  to  be  implied.  Tests  using  this  condition  have  been 
conducted  and  the  results  appear  to  be  satisfactory  for  a  limited  time  after  impingement  of 
the  wave  train  on  the  boundary.  Demonstration  of  this  is  shown  in  a  worked  example  below. 

Spatial  Integration  of  the  Governing  Equations 

The  GN  governing  equations  are  partial  differential  equations,  quasi-linear 
second-order  in  time  and  third-order  in  space.  Further,  these  equations  are  subject  to 
so-called  two-point  boundary  conditions:  a  wave  maker  boundary  condition  at  i  =  2  and 
either  reflection  or  open  boundary  conditions  at  i  =  ns  +  2  or  ns  +  3,  respectively.  It  is 
mathematically  convenient  to  deal  with  the  spatial  problem  first  and,  in  this  regard,  Thomas’ 
algorithm  is  easily  adapted  for  this  purpose  (as  described  in  Demirbilek  and  Webster  (1992)). 
This  algorithm  is  contained  in  subroutine  solve  (see  Figure  2).  Execution  of  the  algorithm 
requires  several  vector  and  matrix  operations  and  these  are  performed  by  separate  simple 
subroutines  as  shown  in  this  figure.  The  coefficients  which  define  the  parti^  differential 
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equations  are  evaluated  in  another  subroutine  called  coefp. 


Temporal  Integration  of  the  Governing  Equations 

The  output  from  solve  is  the  first  time  derivatives  of  the  three  independent  variables 
at  each  point  in  the  computational  domain.  The  value  of  the  variables  at  the  next  time  step  is 
determined  by  using  a  modified  Euler  integration  which  is  second-order  accurate.  The  flow 
near  the  wave  maker  and  near  the  right-hand  boundary  (if  an  open  boundary  is  used)  suffers 
from  an  alternate  point  instability.  That  is,  if  not  corrected,  the  values  of  the  three 
independent  variables  tend  to  oscillate  high  and  low  at  neighboring  points  on  the  finite 
difference  grid.  A  modest  smoothing  of  the  first  few  and  the  last  few  points  in  the 
computational  domain  is  sufficient  to  eliminate  this  cosmetic  blemish.  This  smoothing  is 
accomplished  in  the  subioutiMe  filter.  For  programming  reasons,  this  subroutine  is  called 
before  solve  is  called. 

Finally,  whenever  output  is  performed,  three  additional  dependent  variables,  p,  p., 
and  p,  are  computed  in  the  subroutine  press.  These  variables  are  the  pressure  in  the  fluid 
field  at  the  pressure  gage,  the  integrated  pressure  (from  the  ocean  bottom  to  the  crest  of  the 
wave)  and  the  first  moment  of  this  pressure  (taken  about  the  undisturbed  free  surface  level). 

The  integration  continues  until  either  the  maximum  specified  time,  tmax,  for  the 
simulation  is  reached  or  until  a  wave  breaks  (this  is  exhibited  by  a  failure  of  the  matrix 
inversion  in  the  Thomas  algorithm). 
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PART  TV:  MODEL  INPUT 


The  input  file  communicates  to  the  program  the  type  of  problem  to  be  solved,  the 
wave  system  to  be  created  by  the  wave  malcer,  the  output  desired,  and  the  geometry  of  the 
numerical  wave  tank. 

Before  the  former  quantities  are  discussed,  it  is  necessary  to  specify  the  geometry  of 
the  wave  tank.  The  bottom  of  the  tank  does  not  need  to  be  flat,  but  can  be  a  complex  shape. 
On  the  other  hand,  typical  problems  may  involve  a  spatial  discretization  with  thousands  of 
points  and  it  is  not  practical  to  require  input  of  the  local  depth  for  this  number  of  points.  As 
a  result,  a  simplification  was  adopted.  The  bottom  profile  is  input  as  a  polygonal  shape, 
formed  as  a  series  of  straight  lines  connecting  a  sequence  of  arbitrarily  spaced  nodes.  Figure 
3  shows  this  geometry. 


oomputaSonal  domain 


(xbottm(1).  abottom(1)] 


[xbottm(nbottm),  abottom(nbottm)] 


Figure  3.  Schematic  of  the  nodes  of  the  polygonal  bottom  description 


The  subroutine  bottom  uses  this  polygonal  shape  to  determine  the  water  depth,  a(x), 
at  each  of  the  individual  finite  difference  nodes  (from  i  =  1  to  i  =  ns  +  4).  If  the  value  of 
X  for  a  given  node  falls  outside  of  the  range  of  bottom  data,  then  a(x)  is  assumed  to  be  the 
same  as  that  of  the  closest  input  bottom  data  point.  If  the  value  of  x  falls  within  the 
specified  bottom  data,  a  linear  interpolation  between  the  two  closest  bottom  data  points  is 
used.  Since  the  polygon  does  not  have  even  second  derivatives  at  the  corners  and  since  the 
theory  involves  derivatives  up  to  and  including  the  third  derivative  of  o(x),  the  final  set  of 
interpolated  values  of  o(x)  is  smoothed  using  a  numerical  filter. 
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An  attempt  was  made  to  make  the  input  as  simple  as  possible  without  compromising 
the  ability  of  the  program  to  do  useful  analyses.  The  “input”  file  is  a  flat  text  file.  The 
name  of  this  file  is  not  important,  as  the  user  will  be  ask^  by  the  program  to  specify  it. 
Figure  4,  below,  shows  the  input  for  a  sample  problem  shown  later  in  this  report.  Input 
information  in  this  file  is  confined  to  the  first  60  columns  of  the  file  and  the  remaining 
columns  are  used  for  comments,  a  practice  which  is  highly  recommended.  The  “input”  file 
contains  the  following  information^  : 


Collision  Against  Hydraulic  Gate 
English 


Title 

units  (En  or  SI) 


reflection 

rhs  b.c. 

(op  or  re) 

1200.00000 

x-domain 

180.00000 

t-domain 

13 

#  of  snapshots 

snapl 

60.00000 

snapshot 

title 

& 

time 

snap2 

70.00000 

snapshot 

title 

& 

time 

snap3 

80.00000 

snapshot 

title 

& 

time 

enap4 

90.00000 

snapshot 

title 

& 

time 

snap5 

100.00000 

snapshot 

title 

& 

time 

snap6 

110.00000 

snapshot 

title 

& 

time 

enap7 

120.00000 

snapshot 

title 

& 

time 

snap8 

130.00000 

snapshot 

title 

& 

time 

snap9 

140.00000 

snapshot 

title 

& 

time 

snaplO 

150.00000 

snapshot 

title 

& 

time 

snapll 

160.00000 

snapshot 

title 

& 

time 

8napl2 

170.00000 

snapshot 

title 

6 

time 

snapl3 

180.00000 

snapshot 

title 

& 

time 

3 

#  of 

gages 

gagel  800.0000  -27.00000  0.00000  180.00000 
gage2  1000.0000  -27.00000  0.00000  180.00000 
gage3  1205.0000  -8.00000  0.00000  180.00000 


gage  title,  x,  depth  &  times 
gage  title,  x,  depth  &  times 
gage  title,  x,  depth  &  times 
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#  nodes  in  profile 

0.00000  -30.00000 

profile 

node 

1000.00000  -30.00000 

profile 

node 

1005.00000  -29.93750 

profile 

node 

1010.00000  -29.75000 

profile 

node 

1015.00000  -29.43750 

profile 

node 

1020.00000  -29.00000 

profile 

node 

1180.00000  -11.00000 

profile 

node 

1185.00000  -10.56250 

profile 

node 

1190.00000  -10.25000 

profile 

node 

1195.00000  -10.06250 

profile 

node 

1200.00000  -10.00000 

profile 

node 

1 

#  wave  components 

4.00000  4.00000 

height, 

,  period 

Figure  4.  Sample  input  file 


^  The  internal  FORTRAN  name  is  shown  in  Courier  type  at  the  end  ot  each  description,  and  the 
FORTRAN  format  specification  is  shown  in  parentheses  at  the  end  of  the  last  line  of  each  description. 
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line  1; 


A  title  for  the  computation  which  is  reproduced  as  the  first  line  of  each  output 
file,  title.  (a60) 


line  2:  The  specification  of  the  units:  English  units  &  SI  units  are  supported.  Only 

the  first  two  letters  are  interpreted,  English  or  En  are  equally  acceptable  (there 
is  no  abbreviation  of  SI).  All  subsequent  input  and  output  quantities  are 
assumed  to  be  in  these  units,  qunite.  (a2) 

line  3:  The  specification  of  the  type  of  shore  side  boundary:  open  or  reflection.  Only 

the  first  two  letters  are  interpreted  and  therefore  one  can  write  op  for  open  or 
re  for  reflection,  rhsbdy.  (a2) 

line  4:  The  horizontal  extent  of  the  computational  domain,  xmax.  ( f  lo .  5 ) 

line  5:  The  maximum  time,  tmax  (fio.S) 

line  6:  The  number  of  desired  spatial  snapshots  of  the  variables  (must  be  20  or  fewer 

specifications),  nsnap.  (15) 

Specifications  for  the  snapshots  (one  snapshot  specification  per  line). 

Each  specification  consists  of: 

1.  The  title  of  the  new  file  for  the  output,  tsnap. 

2.  The  time  at  which  the  snapshot  is  to  be  taken,  tss.  (aio,  f lo. 5 ) 

next  line:  The  number  of  desired  wave  gage  time  histories  (must  be  20  or  fewer 

specifications),  ngage.  (15) 

Specifications  for  the  time  histories  (one  time  history  per  line).  Each 
specification  consists  of: 

1.  The  title  of  the  new  file  for  the  output,  tgage. 

2.  The  x-location  of  the  wave  gage,  xgage. 

3.  The  depth  of  the  wave  gage,  gagedp. 

4.  The  starting  time,  tsgi. 

5.  The  ending  time,  t8g2.  (aio,4fio.5) 

next  line:  The  number  of  nodes  in  the  bottom  polygon  (must  be  50  or  fewer  nodes), 

nbottm.  (15) 

next  The  coordinates  of  the  next  node  of  the  bottom  polygon.  The  first  node 

nbottm  must  correspond  to  the  location  of  the  wave  maker;  the  last  node  must 
lines:  correspond  to  the  location  where  an  open  boundary  condition  is  enforced  (if 

used),  (xbottm,  abottm) .  (2fl0.5) 

next  line:  The  number  of  wave  components  for  the  wave  maker,  nwaves.  ( 15 ) 

If  nwaves  ■  1,  then 


next 

ngage 

lines: 


next 

nsnap 

lines: 
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next  line: 


next 

nwaves 

lines: 


The  specifications  of  a  regular  wave: 

1.  wave  height,  height. 

2.  wave  period,  period.  (2fio 

If  nwaves  >  1,  then 

The  components  of  a  random  wave  train  (one  line  per  component): 

1.  wave  amplitude,  ampwav. 

2.  wave  frequency,  freqwv. 

3.  wave  phase.  If  the  phase  of  the  wave  is  larger  than  180',  then  the 
phases  will  be  replaced  by  random  phases,  phaswv.  (3£io.5) 


PARTY;  MODEL  OUTPUT 


There  are  three  types  of  output  files:  an  echo  of  the  original  data  (written  into  the 
"output"  file),  snapshot  files,  and  wave  gage  files. 

The  input,  together  with  any  error  messages  and  a  few  derived  quantities,  such  as  the 
number  of  time  steps  and  the  number  of  spatial  steps,  is  echoed  to  the  computer  screen. 

Since  the  information  on  the  screen  is  volatile,  this  information  (except  for  fatal  error 
messages)  is  also  saved  in  an  output  file  for  future  reference.  The  purpose  of  this  output  file 
is  to  provide  evidence  that  the  input  quantities  were  indeed  read  in  correctly  and  to  form  a 
permanent  record  of  the  calculation. 

The  snapshot  files  have  a  header  describing  the  computation,  followed  by  a  tabular 
list  of  values  of  the  independent  and  dependent  variables,  'hiese  values  are  given  for  the 
whole  computational  domain,  plus  at  the  four  points:  i  =  1,  i  =  2,  i  =  ns  +  3,  and  i  =  ns 
+  4.  These  latter  points  are  outside  the  computational  domain  and  are  needed  for  restart. 
Figure  S  is  a  remnant  of  a  snapshot  file  corresponding  to  the  input  file  given  above.  (The 
numbers  in  this  remnant  were  reduced  by  the  least  significant  figure  so  that  it  can  be  printed 
on  the  page). 

The  header  contains  the  title  on  the  first  line  and  the  second  line  is  blank.  The  third 
line  contains  the  specified  units  and  the  fourth  line  the  specified  right-hand  boundary 
condition.  The  fifth  line  is  blank.  The  sixth  line  contains  the  exact  time  of  the  snapshot 
(which  may  be  slightly  different  from  that  specified  in  the  input  file,  as  a  result  of  die  time 
discretization).  The  seventh  line  contains  the  time  step  index  corresponding  to  the  time  and 
the  eighth  line  contains  dx,  the  spacing  between  nodes.  The  ninth,  tenth,  and  eleventh  lines 
are  blank.  The  twelfth  line  contains  a  list  of  column  tides  and  the  remaining  lines  (beginning 
at  the  thirteenth  line)  are  output  data.  The  first  two  lines  and  last  two  lines  of  output  data 
contain  information  on  the  points  outside  the  computational  domain.  The  column  titles  and 
the  output  data  columns  are  separated  by  tabs  so  tihat  the  entire  matrix  can  be  easily  pasted 
into  a  spread  sheet  or  into  popular  graphing  programs. 

The  meanings  of  the  variables  in  the  output  files  are: 

location: 

X  The  distance  from  the  first  point  in  the  computational  domain, 
independent  variables: 

beta  The  amplitude  of  the  wave  above  the  undisturbed  free  surface. 

uo  The  component  of  velocity  which  is  constant  in  z. 

ui  The  component  of  velocity  that  varies  linearly  with  z. 
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dq)endent  variables: 

pbar  The  pressure  measured  at  the  wave  gage  (at  a  z  coordinate  of  gagedp) . 

po  The  integrated  pressure  (integrated  from  the  bottom  to  the  actual  water 
surface  at  this  location. 

pi  The  first  moment  of  pressure  (the  pressure  multiplied  by  z  and 

integrated  from  the  bottom  to  the  actual  water  surface  at  this 
location. 

All  or  these  variables  are  expressed  in  units  consistent  with  those  given  in  the  input 
specification.  For  instance,  if  English  units  are  selected,  beta  is  in  feet,  uO  in  feet/sec,  ul  in 
1/sec,  pbar  in  Ib/ff,  pi  in  Ib/ft  and  pi  in  ft-lb/ft.  Note  that  the  meaning  of  the  “/ft”  in  pO 
and  pi  is  that  these  are  the  force  and  moment  per  foot  of  wave  tank  in  and  out  of  the  paper 
(i.e.,  in  the  y  or  third  dimension). 

The  wave  gage  files  are  time  histories  of  the  same  variables  displayed  in  the  snapshot 
files.  The  user  can  specify  in  the  input  for  each  gage  the  times  to  start  and  stop  the  recording 
of  the  time  series.  Figure  6  shows  a  remnant  of  a  wave  gage  file. 
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Figure  5.  Typical  snapshot  file 
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Figure  6.  Typical  wave  gage  file 


The  actual  usage  of  the  program  is  quite  simple.  Upon  startup,  the  user  is  prompted 
for  two  pieces  of  information:  ^e  name  of  the  input  file  and  the  name  of  the  output  file  (for 
the  echo  of  the  input  data).  The  program  does  not  assume  a  fixed  name  for  these  two  files, 
so  the  user  can  choose  input  and  output  file  names  that  are  descriptive  of  the  particular  wave 
situation  that  these  files  represent.  All  of  the  rest  of  the  input  information  is  contained  in  the 
referenced  input  file.  The  program  reads  the  data  in  the  input  file,  performs  several 
initialization  steps,  and  then  echoes  the  data  to  the  screen  and  output  file. 

Restart 


After  the  input  file  is  read  the  program  asks:  Is  this  a  restart?  (y/n).  If  the  user  has 
completed  a  calculation  and  wishes  either  to  have  more  snapshots  around  a  critical  event,  to 
have  more  wave  gages  or  to  extend  the  previous  computation,  the  answer  should  be:  v;  if 
this  is  the  first  calculation  with  this  geometry,  the  answer  should  be:  n. 

If  the  answer  is  v,  the  user  will  be  prompted  with:  restart  file  name  =  ?,  The  user  is 
to  input  the  name  of  the  snapshot  file  corresponding  to  the  beginning  time  of  the  new 
computation.  Each  snapshot  file  captures  the  “state”  of  the  flow  and  therefore  contains 
enough  information  for  restart.  The  input  file  elements  that  result  in  output  quantities  can  be 
different  for  the  rerun  from  those  which  were  in  the  original  input  file.  These  are: 


next  ns  nap 
next 
next  ngage 


line  1 ,  the  title  of  the  computation. 

line  5,  the  maximum  time. 

line  6,  the  number  of  desired  spatial  snapshots. 

lines,  the  snapshot  specifications. 

line,  the  number  of  desired  wave  gage  time  histories. 

lines,  the  snapshot  specifications. 


All  other  lines  in  the  input  must  be  the  same  as  the  calculation  that  produced  the  referenced 
snapshot  file  for  the  program  to  run  correctly. 

The  ability  to  rerun  the  computation  is  a  powerful  tool  that  allows  the  user  to  refine 
the  information  he  wishes  to  obtain  from  the  program. 


Run  Time 


The  time  for  execution  of  the  program  is  proportional  to  the  product  of  two  numbers: 
the  number  of  spatial  nodes  and  the  number  of  time  steps.  The  program  displays  these  two 
numbers  on  the  screen  during  the  echo  of  the  “input”  and  records  them  in  the  “output”  file. 
The  time  required  will  depend,  of  course,  on  the  particular  computer  used  to  perform  the 
calculations.  It  is  suggested  that  the  user  perform  a  computation  and  determine  the 
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“calibration  factor”  for  his  own  system. 

Because  some  large  problems  may  take  hours  to  run,  the  iteration  number  is  printed 
out  at  the  beginning  of  each  time  step.  This  assures  the  operator  that  the  calculation  is 
proceeding  normally.  Further,  an  annotation  on  the  screen  is  made  each  time  a  snapshot  file 
is  written. 

Run  time  for  the  model  varies  depending  on  type  of  simulation  performed  (regular 
versus  random  waves),  on  the  desired  size  of  the  modeling  domain,  and  on  the  computer 
platform  used  to  run  the  program.  On  Cray  YM-P  computer,  the  run  time  for  a  pseudo¬ 
random  wave  of  20  components  (i.e.  spectral  wave  simulation)  was  less  than  four  hours  for  a 
region  extending  10  miles  in  on-offshore  direction.  For  a  single  component  wave  (i.e., 
regular  wave),  the  simulation  time  for  a  4  miles  stretch  was  23  minutes  on  a  386/33  MHz  PC 
and  7.4  minutes  on  a  VAX  8800  series  mainframe  computer. 
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APPENDIX  A:  EXAMPLES 


Example  1:  Reflection  of  Waves  from  a  Spillway  Hydraulic  Gate 


A  schematic  (with  the  vertical  scale  exaggerated)  of  the  geometry  of  a  body  of  water 
downstream  of  a  hydraulic  structure  (say,  a  sluice  gate)  is  shown  in  Figure  Al.  The  water 
depth  is  30  ft  downstream  and  a  smooth  transition  is  provided  in  the  2(X)  ft  in  front  of  the 
sluice  gate  to  a  water  depth  of  10  ft  (at  the  sluice  gate).  GNWave  was  run  to  simulate  the 
situation  were  4-ft-high  waves  of  4-sec  period  progress  from  the  deepwater  region  up  to  the 
gate  and  are  reflected  by  the  gate.  Engineering  questions  for  this  situation  include: 

•  How  high  will  the  waves  be  at  the  sluice  gate? 

•  What  forces  will  be  exerted  on  it? 

•  Where  is  the  center  of  pressure  of  these  forces? 

The  computational  domain  selected  is  1,2(X)  ft  in  length  (as  shown  in  Figure  A2). 

This  provides  for  a  large  number  of  waves  in  the  numerical  tank  and  allows  for  a  long  period 
of  simulated  time  to  transpire  before  the  waves  reflected  by  the  sluice  gate  impact  the  wave 
maker.  The  waves  in  the  deep  part  of  the  tank  are  relatively  short  (about  80  ft)  and  therefore 
behave  more  like  deepwater  waves.  These  same  waves  near  the  sluice  gate  (where  the  water 
depth  is  only  10  ft)  behave  like  shallow-water  waves.  The  input  file  is  shown  as  Figure  4  in 
the  main  body  of  this  report. 


Results  and  Discussion 

Figures  A2  through  A4  show  snapshots  of  the  water  surface  at  10-sec  intervals. 

Several  features  are  apparent.  Before  the  wave  train  reaches  the  sluice  gate,  the  leading  wave 
is  noticeably  higher  than  those  which  follow.  This  is  a  direct  result  of  the  dispersive 
character  of  these  finite  water  depth  waves  (Mei  1989).'  When  this  large  wave  impacts  the 
sluice  gate,  it  produces  larger  waves  (over  10  ft  in  height)  there  as  well.  This  wave  is 
reflected,  the  subsequent  waves  in  the  train  result  in  waves  which  average  8  ft  in  height  and 
whose  envelope  appears  to  oscillate  over  a  period  of  about  five  waves.  The  snapshot  at  t  = 
180  clearly  shows  that  the  reflected  waves  are  beginning  to  impinge  on  the  wave  maker. 
Further,  an  examination  of  the  wave  closest  to  the  sluice  gate  shows  the  beginning  of  an 
instability  resulting  from  this  impingement.  Thus,  continuation  of  the  simulation  is  suspect. 
We  should  note  that  a  similar  situation  would  be  encountered  in  a  physical  wave  tank  of 
limited  size. 

Figure  A5  shows  the  time  history  of  the  wave  amplitude  at  the  three  gage  locations  x 
=  800  ft,  X  =  1,000  ft  and  x  =  1,200  ft  (i.e.,  at  the  sluice  gate).  The  effect  of  the  wave 
reflection  is  to  form  standing  waves,  which  can  be  clearly  seen  by  comparison  of  the  three 
gages.  In  particular,  gage  2  is  very  nearly  at  a  node.  Figure  A6  shows  the  time  history  of 
the  pressure  measured  by  pressure  gages  placed  2  ft  above  the  bottom.  This  history 
resembles  the  wave  elevation  time  history,  but  frequency  doubling  can  be  seen  when  standing 
waves  are  formed. 
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See  References  at  the  end  of  the  main  text. 


Figure  A7  shows  the  history  of  the  total  force  per  foot  of  sluice  gate  (the  variable  p.  in 
Demirbilek  and  Webster  (1992))  and  the  location  of  the  center  of  pressure  (the  ratio  p,/p„, 
Demirbilek  and  Webster  (1992)).  Wave  reflection  by  a  sluice  gate  produces  large  standing 
waves  in  the  neighborhood  of  the  gate  and  large  oscillatory  forces  on  the  gate. 

The  effect  of  the  sloping  seabed  at  the  toe  of  the  structure,  which  may  also  be  treated 
as  a  substantial  berm  in  front  of  the  sluice  gate,  can  readily  be  studied  with  the  model 
GNWave.  For  instance,  wave  results  may  be  produced  for  varying  slopes  to  discern  slope 
effects  in  design  estimates,  including  estimating  wave  forces.  It  may  not  be  possible  to 
estimate  such  effects  from  the  results  of  a  single  run. 
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Figure  Al.  Configuration  for  the  test  of  wave  collision  with  a  sluice  gate 
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Figure  A2.  Waves  approaching  sluice  gate  (shown  at  10-sec  intervals) 
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Figure  A4.  Waves  approaching  sluice  gate  (shown  at  10-sec  intervals) 
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Wave  elevation  measured  at  three  gage  locations 


at  X  ■  800' 


Figure  A6.  Pressure  measured  2  ft  off  bottom  at  three  gage  locations 


of  pressure  time  history  at  the  sluice  gate 


Example  U;  Wave  Transformation  over  an  Open-ended  Sloping  Beach 


The  implementation  of  an  open  boundary  condition  (a  condition  that  allows  the  free 
passage  of  waves)  is  a  matter  of  some  experimentation,  since  there  is  no  exact  theoretical 
statement  for  this  condition.  The  program  GNWave  implements  a  Schrodinger  condition  in 
which  the  celerity  of  the  wave  is  assumed  known.  For  most  situations  this  boundary 
condition  is  adequate,  although  it  can  be  expected  to  fail  eventually. 

For  this  test  of  the  open  boundary  condition,  the  short  numerical  wave  tank  is  200  ft 
long  and  a  straight  sloped  bottom  is  chosen  (see  Figures  A8  and  A9).  Waves  of  5 -sec  period 
and  4  ft  in  height  are  generated  at  the  left-hand  side.  These  waves  are  about  120  ft  long  and, 
therefore,  quickly  fill  the  tank.  The  computation  was  allowed  to  proceed  for  50  sec, 
corresponding  to  the  creation  of  10  waves.  A  single  wave  gage  was  placed  at  x  =  100  ft,  the 
center  of  the  tank.  The  input  file  is: 


Test  of  open  boundary  condition  Title 


English 

units  (En  or  SI) 

open 

rhs  b.c.  (op  or  re) 

200.00000 

x-domain 

50.00000 

t-domain 

5 

#  of  snapshots 

snapl 

10.00000 

snapshot  title  &  time 

8nap2 

20.00000 

snapshot  title  &  time 

snapa 

30.00000 

snapshot  title  &  time 

8nap4 

40.00000 

snapshot  title  &  time 

snaps 

50.00000 

snapshot  title  &  time 

1 

#  of  gages 

gagel  100.00000  -30.00000 

0.00000  50.00000  gage  title,  x,  depth  fi  times 

2 

#  nodes  in  profile 

0.00000 

-30.00000 

profile  node 

200.00000 

-15.00000 

profile  node 

1 

#  wave  components 

4.00000 

5 . 00000 

height, period 

Figure  A8.  Input  for  test  of  the  open-boundary  condition 


Results  and  Discussion 

The  time  history  of  the  surface  elevation  at  the  wave  gage  is  shown  in  Figure  A 10. 

The  first  few  waves  passing  the  gage  were  somewhat  different  in  amplitude  from  the  waves 
which  are  ultimately  achieved  in  steady  state.  Figure  All  shows  five  snapshots  at  intervals  of 
10  sec.  Since  the  wave  period  is  5  sec,  these  snapshots  should  be  identical  when  steady  sjite 
is  achieved.  The  first  snapshot  is  considerably  dif^ferent  from  the  others  because  the  first 
wave  has  not  made  it  to  the  boundary. 

The  differences  between  the  remaining  snapshots  are  more  difficult  to  discern,  but  the 
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differences  can  be  most  clearly  seen  at  the  right-hand  edge  of  the  domain.  The  difference 
between  the  snapshot  at  40  sec  and  that  at  SO  sec  is  less  than  a  line  width.  These  waves  have 
sharper  crests  than  troughs  due  to  nonlinear  effects. 

The  shortening  of  the  wave  in  the  shallower  water  at  the  right-hand  side  of  the  domain 
is  also  obvious  in  these  snapshots.  The  open  boundary  condition  appears  to  work 
satisfactorily  and  permits  over  10  waves  of  moderate  steepness  to  pass  through  the  boundary 
with  no  visible  reflections  or  distortions. 
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Figure  A9.  Configuration  for  test  of  open  boundary  condition 
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When  finite  waves  pass  over  a  submerged  disturbance,  such  as  a  sandbar,  they 
steepen.  If  the  sandbar  is  not  too  pronounced,  these  steep  waves  will  not  break.  When  they 
pass  the  sandbar,  the  waves  once  again  decrease  in  amplitude.  An  interesting  question  is:  are 
these  resulting  waves  substantially  different  in  character  from  the  waves  before  the  sandbar? 
If  the  waves  were  small  enough  so  that  linear  wave  theory  applies  (both  in  the  neighborhood 
of  the  sandbar  and  elsewhere)  then  the  answer  is  clearly  that  there  is  no  remnant  of  the  sand¬ 
bar  except  for  a  change  of  phase  of  the  wave  corresponding  to  the  slower  phase  speed  of  the 
wave  in  the  neighborhood  of  the  bar.  In  the  case  of  waves  steep  enough  to  exhibit 
nonlinearities,  the  authors  know  of  no  literature,  experimental  or  otherwise,  which  describes 
the  situation.  The  case  of  interaction  of  multiple  wave  components  in  this  situation  is  even 
more  complicated. 

In  this  worked  example  a  pseudo-random  wave  train  is  allowed  to  evolve  over  two 
different  bottom  topographies:  the  first  is  water  of  uniform  depth  (in  this  case,  a  10-ft  depth), 
and  the  second  is  water  that  is  10  ft  in  depth  everywhere  except  for  a  “sandbar.”  Each 
numerical  wave  tank  is  800  ft  long  and  an  open  boundary  condition  is  applied  on  the 
right-hand  boundary.  The  sandbar  begins  at  a  distance  of  425  ft  from  the  wave  maker  and 
extends  for  150  ft  downstream.  The  sandbar  is  therefore  centered  about  a  point  500  ft  from 
the  wave  maker.  The  shape  of  the  sandbar  is  a  trapezoid  with  a  flat  top  50  ft  in  length.  The 
pseudo  random  wave  train  is  composed  of  five  waves  each,  with  an  amplitude  of  0.5  ft  and 
each  with  a  different  frequency.  Sketches  of  these  two  situations  are  shown  in  each  of  the 
following  graphs. 

Wave  gages  for  each  situation  are  placed  at  250,  5(X),  and  750  ft  from  the  wave 
maker.  For  the  case  of  the  sandbar,  these  gages  correspond  to  a  point  well  upstream  of  the 
sandbar,  right  over  the  middle  of  the  sandbar,  and  well  downstream  of  the  sandbar.  The  two 
input  files  are  given  below  in  Figures  A12  and  A13. 


Test  of  Random  Waves  in  a  Channel 

Title 

English 

open 

800.00000 

135.00000 

14 

snapO 

30.00000 

units  (En  or  SI) 

rhs  b.c.  (op  or  re) 

x-domain 

t-domain 

#  of  snapshots 

snapshot  title  &  time 

snapl 

40.00000 

snapshot  title  &  time 

Bnap2 

80.00000 

snapshot  title  &  time 

snap2a 

85.00000 

snapshot  title  &  time 

snap2b 

90.00000 

snapshot  title  &  time 

snap2c 

95.00000 

snapshot  title  &  time 

snap2d 

100.00000 

snapshot  title  6  time 

snap2e 

105.00000 

snapshot  title  &  time 

snap2  f 

110.00000 

snapshot  title  &  time 

snap2g 

115.00000 

snapshot  title  &  ti'ne 

snap3 

120.00000 

snapshot  title  &  cime 

snap3a 

125.00000 

snapshot  title  &  time 

snap3b 

130.00000 

snapshot  title  &  time 

snap3c 

135.00000 

snapshot  title  &  time 

3 

gagel  250 

.00000  -10.00000  0.00000 

200.00000 

#  of  gages 

gage  title,  x,  depth  &  times 

gage2  500 

.00000  -10.00000  0.00000 

200.00000 

gage  title,  x,  depth  &  times 

gage3  750 

.00000  -10.00000  0.00000 

200.00000 

gage  title,  x,  depth  &  times 

2 

0.00000 

-10.00000 

/  nodes  in  profile 
profile  node 

800.00000 

-10.00000 

profile  node 

5 

0.50000 

1.00000  200.00000 

#  wave  components 
amplitude,  freq  6  phase 

0.50000 

1.21000  200.00000 

amplitude,  freq  6  phase 

0.50000 

1.33000  200.00000 

amplitude,  freq  6  phase 

0.50000 

1.41000  200.00000 

eunplitude,  freq  &  phase 

0.50000 

1.59000  200.00000 

amplitude,  freq  6  phase 

Figure  A 12.  Input  for  pseudo-random  waves  in  a  uniform  depth  channel 
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Test  of  Random  Waves  over  a  Sandbar 

English 

open 

800.00000 

135.00000 


14 

snapO 

30.00000 

snapl 

40.00000 

snap2 

80.00000 

snap2a 

85.00000 

snap2b 

90.00000 

snap2c 

95.00000 

snap2d 

100.00000 

snap2e 

105.00000 

snap2f 

110.00000 

snap2g 

115.00000 

snap3 

120.00000 

snap3a 

125.00000 

snap3b 

130.00000 

snap3c 

O 

135.00000 

w 

gagel 

250.00000  -10.00000 

0.00000 

gage2 

500.00000  -  6.00000 

0.00000 

gage3 

750.00000  -10.00000 

0.00000 

6 


0.00000  -10.00000 
425.00000  -10.00000 
475.00000  -  6.00000 
525.00000  -  6.00000 
575.00000  -10.00000 
800.00000  -10.00000 
5 

0.50000  1.00000  200.00000 
0.50000  1.21000  200.00000 
0.50000  1.33000  200.00000 

0.50000  1.41000  200.00000 

0.50000  1.59000  200.00000 


Title 

units  (En  or  SI) 
rhs  b.c.  (op  or  re) 
x-domain 
t-domain 
/  of  snapshots 
snapshot  title  &  time 
snapshot  title  &  time 
snapshot  title  £  time 
snapshot  title  fi  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  time 
snapshot  title  £  tiise 

#  of  gages 

200.000  gage  title,  x,  depth  Ji  times 

200.000  gage  title,  x,  depth  £  times 

200.000  gage  title,  x,  depth  £  times 

#  nodes  in  profile 
profile  node 
profile  node 
profile  node 
profile  node 
profile  node 
profile  node 

#  wave  components 
amplitude,  freq  £  phase 
amplitude,  freq  £  phase 
amplitude,  freq  £  phase 
amplitude,  freq  £  phase 
amplitude,  freq  £  phase 


Figure  A 13.  Input  for  pseudo-random  waves  over  a  sandbar 
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Results  and  Discussion 


Two  problems  of  practical  interest  are  considered  in  this  example  and  simulation 
results  for  both  are  jointly  displayed  and  compared  to  illustrate  the  effects  of  a  sandbar  on  the 
down-wave  propagating  past  the  seabed  mound. 

Figures  A 13  through  A19  show  the  profiles  of  the  wave  trains  at  various  times  in  the 
interval  from  0  -  135  sec.  These  figures  are  arranged  so  that  the  top  two  graphs  are  the 
waves  in  the  channel  and  the  bottom  two  are  waves  over  a  sandbar  (a  little  sketch  is  placed  on 
the  graphs  to  differentiate  the  two).  From  these  graphs,  it  is  clear  that  the  sandbar  has  a 
dramatic  effect,  especially  on  the  waves  that  were  initially  steeper.  It  is  also  clear  that 
whereas  the  waves  upstream  of  the  sandbar  are  nearly  identical,  those  downstream  are  quite 
different  in  character. 

Figures  A20  through  A22  show  the  time  histories  of  the  three  wave  gages  for  the  two 
situations.  It  can  be  clearly  seen  that  the  waves  upstream  of  the  sandbar  are  identical.  Those 
over  the  sandbar  demonstrate  the  sharp,  peaked  character  of  waves  in  shoaling  situations. 

Most  interesting  are  the  waves  downstream  of  the  sandbar  (Figure  A21).  These  waves  are 
clearly  much  different  in  nature  from  the  corresponding  waves  in  the  uniform  depth  channel. 
The  peaked  character  of  the  near  breaking  waves  near  the  sandbar  is  retained,  resulting  in  a 
wave  train  which  has  more  high-frequency  content  than  those  in  the  channel. 
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Figure  A14.  Snapshots  of  wave  profiles  at  t  =  30  sec  and  t  =  40  sec 
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Figure  A 16.  Snapshots  of  wave  profiles  at  t  =  90  sec  and  t  =  95  sec 
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Figure  A 17.  Snapshots  of  wave  profiles  at  t  «  100  sec  and  t  =  105  sec 
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Figure  A18.  Snapshots  of  wave  profiles  at  t  =  110  sec  and  t  =  115  sec 
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Figure  A19.  Snapshots  of  wave  profiles  at  t  =  120  sec  and  t  =  125  sec 
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Figure  A23.  Time  histories  of  wave  elevation  at  x  =  750  ft 


Example  IV:  Shoaling  Waves  over  a  Uniformly  Sloped  Coastal  Beach 


As  waves  approach  a  beach,  they  steepen,  become  unsymmetrical  and  eventually 
break.  The  GN  Level  11  theory  which  forms  the  basis  of  GNWave  is  not  capable  of 
predicting  the  evolution  of  waves  beyond  the  breaking  stage,  but  it  can  predict  the  evolution 
almost  up  to  breaking.  The  example  below  is  identical  to  that  presented  in  Shields’  thesis  and 
models  experiments  of  Hansen  and  Svendsen  (1979)  conducted  at  the  Technical  University  of 
Denmark  Regular  waves  are  created  in  a  tank  25  m  long  with  a  beach  that  begins  13  m 
downstream  of  the  wave  maker.  The  beach  has  a  slope  of  1:34.26.  The  water  depth  at  the 
wave  maker  is  0.360  m.  For  the  purposes  of  modeling  here,  the  bottom  of  the  beach  is  not 
allowed  to  pierce  the  water  surface,  but  is  kept  at  a  water  depth  of  0.040  m  from  a  point 
0.02396  m  downstream  of  the  beach.  The  input  file  for  this  example  problem  is  listed  below 
in  Figure  A24. 


Hansen  &  Svendsen  Test  #031041 
SI 

open 
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20.00000 
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snapl2.0 

12.00000 

snapl4.0 
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17.50000 

snapl7.6 
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0.00000 

gagelGO 

19.8500  -  0.16000 

0.00000 

gagel20 

•1 

21.2200  -  0.12000 
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#  wave  components 
height, period 


Figure  A24.  Input  file  for  the  shoaling  of  regular  waves 
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A  small  sketch  of  the  test  situation  is  shown  at  the  top  of  Figure  A25.  Regular  waves  are 
generated  at  the  wave  maker  with  a  height  of  40  mm  and  a  period  of  3.33  sec  (frequency  of 
0.3  Hz). 


Results  and  Discussion 

Figure  A25  shows  five  wave  profiles  at  times  from  17.2  sec  to  18.4  sec  (the  first  wave 
breaks  at  18.5  sec).  Figure  A26  shows  the  detail  of  the  waves  near  the  beach  at  0.1-sec 
intervals.  For  these  details,  the  output  data  (shown  as  small  squares)  were  fitted  with  spline 
curves  to  yield  smooth  wave  shapes.  The  sharpening  and  steepening  of  the  wave  is  apparent, 
as  is  the  development  of  a  “shoulder  wave”  just  behind  the  main  crest. 

These  results,  although  similar  in  form  to  those  presented  by  Shields  (1986)  are 
different  in  some  respects.  In  particular,  the  waves  break  further  up  the  beach  from  the  point 
reported  by  Shields  and  that  measured  by  Hansen  and  Svendsen  (1979).  An  investigation  of 
this  was  made  and  it  was  determined  that  these  discrepancies  were  due  to  the  built-in  facilities 
in  GNWave  which  yield  an  easy  user  interface  and  do  not  reflect  on  either  the  theory  or  its 
general  usage. 

The  first  wave  in  the  train  is  purposely  attenuated  (or  “ramped”)  and  this  is 
appropriate  for  most  problems.  In  this  problem,  the  first  full  wave  breaks  and  the 
computation  cannot  continue.  Thus,  although  the  subsequent  waves  in  the  wave  train  have  the 
required  height  of  40  mm,  the  first  wave  is  not  quite  this  high.  This  causes  the  wave  to  break 
further  up  the  beach  from  the  point  determined  by  Shields  and  that  measured  by  Hansen  and 
Svendsen  (1979).  Shields  (1986)  was  able  to  get  around  this  problem  by  starting  the 
calculation  with  an  already  developed  train  of  regular  waves  and  letting  them  progress  on 
shore.  This  way  his  first  wave  was  at  the  full  height  and  he  was  able  to  compare  directly 
with  the  experimental  results.  It  is  not  possible  to  allow  such  tinkering  with  the  input  waves 
and  still  have  an  easy-to-use  program. 

Further  validation  of  GNWave  with  field  and  laboratory  data  is  in  progress. 
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Figure  A26.  Wave  profiles 
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